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Abstract 


A conical Euler code was developed to study unsteady vortex- dominated 
flows about rolling, highly swept delta wings undergoing either forced 
motions or free-to-roll motions that include active roll suppression. The 
flow solver of the code involves a multistage, Runge-Kutta time-stepping 
scheme that uses a cell- centered, finite-volume, spatial discretization of 
the Euler equations on an unstructured grid of triangles. The code 
allows for the additional analysis of the free-to-roll case by simultaneously 
integrating in time the rigid-body equation of motion with the governing 
flow equations. Results are presented for a delta wing with a 75° swept, 
sharp leading edge at a free- stream Mach number of 1.2 and at 10°, 20° , 
and 30° angle of attack a. At the lower angles of attack (10° and 20°), 
forced-harmonic analyses indicate that the rolling-moment coefficients 
provide a positive damping, which is verified by free-to-roll calculations. 
In contrast, at the higher angle of attack (3(F), a forced-harmonic analysis 
indicates that the rolling-moment coefficient provides negative damping at 
the small roll amplitudes. A free-to-roll calculation for this case produces 
an initially divergent response, but as the amplitude of motion grows 
with time, the response transitions to a wing-rock type of limit cycle 
oscillation, which is characteristic of highly swept delta wings. This limit 
cycle oscillation may be actively suppressed through the use of a rate- 
feedback control law and antisymmetrically deflected leading- edge flaps. 
Descriptions of the conical Euler flow solver and the free-to-roll analysis 
are included in this report. Results are presented that demonstrate how the 
systematic analysis of the forced response of the delta wing can be used 
to predict the stable, neutrally stable, and unstable free response of the 
delta wing. These results also give insight into the flow physics associated 
with unsteady vortical flows about delta wings undergoing forced motions 
and free-to-roll motions, including the active suppression of the wing-rock 
type phenomenon. The conical Euler methodology developed is directly 
extendable to three-dimensional calculations. 


Introduction 

In recent years, the understanding and prediction 
of the complex flows about high-performance aircraft 
at high angles of attack have generated much interest 
within the fluid dynamics community. (See refs. 1 
and 2.) These aircraft typically have thin, highly 
swept lifting surfaces that produce vortical flow over 
the leeward side of the vehicle at high angles of at- 
tack. This vortical flow can have beneflcial effects, 
such as lift augmentation at high angles of attack, 
on performance. However, it also may have adverse 
effects, such as structural fatigue due to tail buffet 
and stability and control problems due to wing rock, 
wing drop, nose slice, and pitch-up. (See ref. 3.) Con- 
sequently, considerable experimental work has been 
done to understand the basic flow physics of vorti- 
cal flows about delta wings at high angles of attack. 


Experimental research efforts directed toward under- 
standing and documenting steady vortical flows are 
typifled by the detailed flow-fleld measurements 
about simple-delta, cranked-delta, and canard-delta 
wing conflgurations at low speed (ref. 4) and the 
low-speed tests on a 75° swept delta wing (ref. 5). 
For supersonic free-stream Mach numbers, vortical 
flows have been measured by Squire (ref. 6) for an 
elliptic cone delta wing and by Miller and Wood 
(ref. 7) for a series of swept sharp- leading-edge delta 
wings. Experiments to investigate unsteady vorti- 
cal flows for forced harmonic and free-to-roll mo- 
tions of an 80° swept delta wing at low speeds have 
been reported in reference 8. This wing underwent 
self-induced periodic roll oscillations known as “wing 
rock” for angles of attack greater than 25°. Levin and 
Katz (ref. 9) tested 76° swept and 80° swept delta 
wings and found that only the 80° model exhibited 




wing rock at high angles of attack. Further studies 
have been performed by Jun and Nelson (ref. 10) and 
Arena and Nelson (ref. 11). These studies show, for 
example, the time histories of the vortex core posi- 
tion during a cycle of wing rock (ref. 10) and the 
static and dynamic effects due to vortex breakdown 
(ref. 11). Also, reference 12 contains experimental 
water-tunnel results that show wing rock for several 
delta wing planforms along with detailed flow visu- 
alization diagrams. Although much work remains to 
be done, references 4 to 12 have contributed signifi- 
cantly to the understanding of steady and unsteady 
vortex-dominated flow fields. 

From a computational point of view, consider- 
able effort also has been spent on developing methods 
for predicting steady and unsteady vortex-dominated 
flows. (See refs. 13 and 14.) Hoeijmakers (ref. 13) 
gives a review of computational methods for the 
determination of steady vortical flow characteristics 
with an emphasis on classical methods such as dis- 
crete vortex, cloud in cell, panel, vortex layer with 
finite core, leading-edge suction analogy, and vortex 
lattice. With respect to unsteady methods, a non- 
linear mathematical model is presented in refer- 
ence 15 for calculating wing-rock characteristics 
based on aerodynamic derivatives that are evaluated 
using steady- flow aerodynamics. Wing rock was sim- 
ulated in references 16 to 18 by using an unsteady 
vortex-lattice method to predict the aerodynamic 
loads, and the equation of rolling motion was inte- 
grated by using a predictor-corrector method. The 
methods of references 15 to 18 were used to predict, 
with reasonable accuracy, the low-speed wing-rock 
characteristics of the delta wings studied in refer- 
ences 8 and 9. Use of the more modern compu- 
tational fluid dynamics techniques for the predic- 
tion of vortex- dominated flows (ref. 14) has focused 
primarily on steady applications (refs. 19 to 26); 
there are notable exceptions where applications have 
been made to rolling delta wings that were under- 
going forced harmonic (refs. 27 to 29) and free-to-roll 
(ref. 30) motions. Kandil and Chuang, for example, 
calculated flows past rolling delta wings by using the 
conical Euler equations for sharp-leading-edge wings 
(ref. 27) and the conical Navier-Stokes equations for 
rounded- leading-edge wings (ref. 28). In reference 29, 
results for a rolling delta wing were computed with a 
conical Euler flow solver on an unstructured grid of 
triangles. The methods of reference 29 were extended 
in references 30 to 32 to include a free-to-roll capa- 
bility, and results are shown for a freely rolling delta 
wing that exhibited a limit cycle or wing-rock type 
motion that is characteristic of highly swept delta 
wings. Subsequent delta- wing calculations obtained 


by using the conical Euler equations on a structured 
mesh also exhibited limit cycle oscillations at high 
angles of attack. (See ref. 33.) 

The objective of the current research is to study 
unsteady, vortex-dominated flow fields by using the 
conical Euler equations as a first step in investigat- 
ing the three-dimensional problem. The purpose of 
this paper is to report on the development of a coni- 
cal Euler analysis method to study unsteady, vortex- 
dominated flows about rolling delta wings under- 
going either pulsed motion, forced harmonic motion, 
or free-to-roll motion that includes active roll sup- 
pression. Descriptions of the conical Euler flow solver 
and free-to-roll analysis are included. The flow solver 
involves a multistage, Runge-Kutta time-stepping 
scheme and a cell-centered, finite- volume, spatial dis- 
cretization of the Euler equations on an unstructured 
grid of triangles. The code was modified to include 
the simultaneous time integration of the rigid-body 
equation of motion with the governing flow equa- 
tions to allow for the additional analysis of the free- 
to-roll case. The analysis also includes a capabil- 
ity for implementing an active feedback control law 
with antisymmetrically deflected leading-edge flaps 
for suppression of the wing-rock motion. Limited ex- 
perimental and numerical work has been conducted 
on the use of flaps (refs. 34 and 35) and leading-edge 
blowing (refs. 36 and 37) for roll control. Results 
are presented herein for a 75° swept, sharp-leading- 
edge delta wing at a free-stream Mach number of 1.2 
and at 10°, 20°, and 30° angle of attack. These 
results demonstrate how the systematic analysis of 
the forced response of the delta wing can be used 
to predict the stable, unstable, and neutrally sta- 
ble free response of the delta wing. These results 
also give insight into the flow physics associated with 
unsteady vortical flows about delta wings that are 
undergoing forced motions and free-to-roll motions, 
including the active suppression of the wing-rock 
type phenomenon. 

Symbols 


aoo free-stream speed of sound 

poo 

C2 

Cl rolling- moment coefficient 

transfer function of rolling- moment 
coefficient due to flap deflection 

transfer function of rolling- moment 
coefficient due to roll 
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Px 

p 

poo 


root chord of wing 
total energy 

mass moment of inertia about 
longitudinal axis 

imaginary unit, y/—l 

control gain 

reduced frequency of oscillation 
based on half the root chord, 

rolling moment, positive clockwise 
when viewed from aft 

free- stream Mach number 

fluid pressure 

total pressure loss 

free-stream dynamic pressure 

planform area 

transpose operator 

dimensional time 

nondimensional time, 

nondimensional time at center of 
pulse 

contravariant velocities in x, y, and 
z directions 

free-stream velocity 

components of velocity vector in x, 
?/, and z directions 

Cartesian coordinates 

angle of attack 

ratio of specific heats 

energy exchange during harmonic 
cycle 

normalized energy exchange during 
harmonic cycle 

dimensionalized global time step 

nondimensional time step 

leading-edge flap deflection angle, 
positive clockwise when viewed from 
aft 

structural damping 
density 

free-stream density 


(f) instantaneous roll angle, positive 

clockwise when viewed from aft 

4>o harmonic and pulse roll- angle 

amplitude 

uj frequency, radians per second 

Primes with symbols indicate differentiation with 
respect to nondimensional time. 


Governing Equations 


The flow is governed by the time-dependent Euler 
equations, which may be written in conservation-law 
form as 


dQ dE d¥ dG 
dt dx ^ dy ^ dz 


( 1 ) 


where Q is the vector of conserved variables defined 
by 

Q = [p, pu, pv, pw, ef (2) 


and E, F, and G are the convective or inviscid fluxes 
given by 


pu 

pUu -h p 

E = I PU'^ 

pUw 

{e-\-p)U -h xtp 

pV 
pVu 

F = j P^'^ + P 
pVw 

(e + p) y + ytp 

pW 
pW u 

G = I pWv 

pWw + p 
,{e-\-p)W + ztp. 


(3a) 


(3b) 


(3c) 


The contravariant velocities C/, V, and W are defined 
by 

U = u — xt\ 


V = v-yt > 


(4) 


W = w ~ zt j 

where x^, yt^ and zt are the grid speeds in the x, 
y, and z directions, respectively. The pressure p is 
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determined by the equation of state for a perfect gas 
as follows: 


p = (7 - 1) 







( 5 ) 


Equations (1) to (5) have been nondimensionalized 
by the free-stream density and the free-stream speed 
of sound. 

If interest is restricted to supersonic flow past con- 
ical bodies, the conical-flow assumption that is exact 
for steady inviscid supersonic flow can be made. This 
assumption reduces the problem from three dimen- 
sions to two dimensions, which significantly decreases 
the computational resources that are required to in- 
vestigate such flows. For unsteady flows, however, 
the conical assumption implies instantaneous prop- 
agation along radial lines. The following change of 
variables is then required: 





(6) 


The three-dimensional Euler equations then reduce 
to 


^ + ^(F-pE) + ^(G-CE) + 2E = 0 (7) 

Equation (7) may be rewritten in integral form for 
solution as 

^ Qdr,d<:+ I |(F - 7?E) dC - (G - <E) dr,] + /2E d7? dC = 0 

Jq JdQ Jn 

( 8 ) 

where is the integration area and the second 
integral is a boundary integral that results from 
application of the divergence theorem. 

Solution Algorithm 

In this section, an algorithm for the solution of the 
unsteady conical Euler equations on an unstructured 
mesh of triangles is described. Conceptually, this 
spatial discretization reduces to central differencing 
on a rectangular mesh. 

Spatial Discretization 

The conical Euler equations in integral form 
(eq. (8)) are solved numerically by using a finite- 
volume algorithm developed for analysis with an un- 
structured grid that is made up of triangles. The 


algorithm is a cell-centered scheme in which the flow 
variables are stored at the centroid of each triangle; 
the control volumes are the triangular elements of 
the mesh. The boundary integral in equation (8) 
is formed by the fluxes E, F, and G evaluated at 
the midpoint of each edge of a given triangular el- 
ement. The inviscid fluxes at the midpoint at each 
edge are evaluated with a cell-centered type of ap- 
proach in which the fluxes of two adjacent elements 
are averaged across a common edge. For example, 
along element edge m-fc in figure 1, the contribution 
to the boundary integral in equation (8) for element 
i is given by 

“ Vmk^mk) (Cm “ Cfc) ~ ~ Cmk^mk) (Vm “ Vk) 

( 9 ) 


(10a) 

(lOb) 

(10c) 

(lOd) 

(lOe) 


where 

Vmk ~ 2 ^k) 

Cmk ~ 2 

^mk = \ (Ei +Ej) 
^mk ~ 2 ^3 ) 

^mk = 2 


Artificial Dissipation 

Artificial dissipation is added explicitly to pre- 
vent oscillations near shock waves and to damp high- 
frequency, uncoupled error modes. Specifically, an 
adaptive blend of harmonic and biharmonic opera- 
tors corresponding to second- and fourth-difference 
dissipation, respectively, is used; this blend is similar 
to the dissipation described in reference 38. The bi- 
harmonic operator provides a background dissipation 
to damp high-frequency errors, and the harmonic op- 
erator prevents oscillations near shock waves. The 
harmonic operator is multiplied by a pressure switch 
that is first-order accurate near shocks and second- 
order accurate in smooth regions of the flow. The 
biharmonic operator is third-order accurate and is 
adaptively turned off to prevent overshoots in regions 
of shock waves. The harmonic-difference operator 
stencil for element i with surrounding neighbors 1, 
2, and 3 is shown in figure 2. For elements adjacent 
to a boundary, image or “ghost” cells are used to 
complete the stencil shown in figure 2. The value? 
of the conserved variables associated with the image 
cells are determined by the type of boundary that the 
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shared edge represents. If the boundary is a solid 
wall, then the normal and tangential velocities are 
set in the image cell so that a free-to-slip condition 
is imposed on the boundary edge. Specifically, the 
normal velocity in the image cell is set opposite to 
that of the adjacent interior cell; tangential velocity 
is simply mirrored, so that when the velocities are 
averaged across the boundary edge, a flow tangency 
condition is enforced. The pressure and density are 
refiected into the image cells such that a zero nor- 
mal derivative is imposed on the boundary edge. If a 
cell is adjacent to a far-field boundary, then the con- 
served variables in the adjacent image cell are set to 
free-stream conditions. This boundary condition re- 
quires that the bow shock be captured within the in- 
terior of the computational domain. The biharmonic 
operator for an element is computed with second- 
difference terms from the harmonic operators. Along 
the boundaries, the second-difference term in the im- 
age cell is simply set equal to its value in the parent 
cell, so that no additional ghost cells are required at 
the boundaries. 

Boundary Conditions 

The boundary conditions are enforced by speci- 
fying the fluxes along boundary edges for use in the 
evaluation of the boundary integral in equation (9). 
For those edges along the far-field boundary, free- 
stream conditions are applied with the assumption 
that the bow shock is captured. To impose a no- flow- 
through boundary along the surface of the body, the 
boundary integral in equation (8) is first rewritten in 
terms of the flux velocity defined by 

{V - 7]U) AC - (VF - CU) Ar/ (11) 

For edges along the solid boundary, this flux velocity, 
which is proportional to the velocity normal to the 
edge, is then set equal to zero. The pressure terms 
along the solid boundary edge are evaluated in the 
boundary integral with the cell-center values. 

Time Integration 

The conical Euler equations are integrated in time 
by assuming that the conserved variables represented 
by Q are constant within a control volume which 
yields 


j^{A,Qi) + C{Qi)-D{Qi) = 0 ( 12 ) 

where C and D are the convective and dissipative 
operators, respectively, and Ai is the area surround- 
ing element i. These equations are integrated in time 


by using an explicit four-stage, Runge-Kutta, time- 
stepping scheme given by 


q(i) 

q ( 2 ) 

q ( 3 ) 

q ( 4 ) 


qn 
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^n+1 

^n+1 

A^ 

q ( 4 ) 


3(0) _ i 


1 At 


(0) 


Q 


q(0), 


3 > 1^+1 
1 At 


2 yl^+1 


[ c ( q «»)- z ?( q ( o ))] 

[c(Qb))-D(Q(«))] 

[c(q(2))-d(q(o))] 


> ( 13 ) 


At 


[C ( q ( 3 >) -£>( q ( 0 ))] 


where At is the global time step and the superscript 
n represents the value at the time level n. In this 
scheme, the convective operator is evaluated at each 
stage; for computational efficiency, the dissipative 
operator is evaluated only at the first stage. The 
Runge-Kutta scheme represented by equations (13) 
is second-order accurate in time and includes the 
necessary terms to account for changes in cell areas 
as the result of a moving or deforming'mesh. 


Implicit Residual Smoothing 

The explicit time-integration scheme described in 
the preceding section has a step size that is lim- 
ited by the Courant-Friedrichs-Lewy (CFL) condi- 
tion corresponding to a CFL number of 2\/2. To 
accelerate convergence to the steady state, the CFL 
number may be increased by averaging the residual 
with values at neighboring elements. This resid- 
ual averaging is accomplished by replacing R^ with 
the smoothed residual R^ given by 

R^ — eV^Rj = R^ (14) 

where c is a constant that controls the amount of 
smoothing and is an undivided Laplacian oper- 
ator. These implicit equations are solved approxi- 
mately by using several Jacobi iterations similar to 
those in references 38 and 39. 

For steady-state calculations, convergence is fur- 
ther accelerated by using enthalpy damping (ref. 40) 
and local time stepping. The local time stepping uses 
the maximum allowable step size at each grid point 
as determined by a local stability analysis. For un- 
steady applications, however, a global time step must 
be used because of the time- accuracy requirement. 
By using a time-accurate version of equation (14) 
similar to that of reference 41, the maximum allow- 
able global time step may be increased to a value 
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that is larger than that dictated by the CFL condi- 
tion. In this procedure, the constant e is replaced by 
a parameter defined by 


€ = max 


1 

4 



0.0 


(15) 


which varies from grid point to grid point. In equa- 
tion (15), At is the time step taken and A^qpl is the 
locally allowable time step for the four-stage Runge- 
Kutta time-stepping scheme. 


Deforming-Mesh Algorithm 

The deforming-mesh algorithm, as developed in 
reference 42, models the triangular mesh as a spring 
network in which each edge of a triangle represents 
a spring with stiffness inversely proportional to the 
square of its length. In this method, the grid points 
along the outer boundary are held fixed, and the 
grid points along the wing (inner boundary) are 
specified. The locations of the interior points are 
then determined by solving the static equilibrium 
equations, which result from a summation of forces 
at each node in the t] and ^ directions. The solution 
of the equilibrium equations is carried out by using 
a predictor-corrector method that first predicts the 
new locations of the interior points by extrapolation 
from the previous time levels and then corrects these 
locations by using several Jacobi iterations of the 
static equilibrium equations. The predictor-corrector 
procedure is relatively efficient, since it requires only 
a few Jacobi iterations to adequately move the mesh. 


Pulse Transfer-Function Analysis 

Generally, unsteady load coefficients can be ob- 
tained by calculating several cycles of a forced- 
harmonic oscillation and using the last cycle of os- 
cillation to determine the load. This process requires 
one' flow-field calculation for each value of reduced 
frequency of interest. In contrast, the unsteady load 
coefficients may be determined for a wide range of 
reduced frequency in a single flow-field calculation 
by using the pulse transfer- function analysis. The 
pulse transfer-function analysis has been used pre- 
viously to determine the generalized aerodynamic 
forces (GAF), which are used in aeroelastic analyses. 
(See refs. 43 and 44.) In the pulse analysis, the un- 
steady load coefficient is computed indirectly from 
the response of the flow field as a result of motion 
that is represented by a smoothly varying, exponen- 
tially shaped pulse. Results computed by using the 
pulse analysis for a pitching flat-plate airfoil were in 
good agreement with parallel linear-theory calcula- 
tions. (See ref. 44.) Reference 44 also shows that 


the GAF airfoils at transonic speeds that were com- 
puted from a pulse analysis were in good agreement 
with the GAF values that were computed by using 
the harmonic method, which tends to verify that the 
analysis is valid for predicting the small perturbation 
response about a nonlinear flow field. These calcula- 
tions verify the accuracy of the pulse analysis. There- 
fore, because of the computational efficiency of the 
pulse transfer-function analysis, the capability was 
implemented within the conical Euler code to calcu- 
late the rolling-moment coefficient due to roll of 
a delta wing. The pulse in roll angle is expressed as 


<t> (t) = <t >0 exp 


-M^ {t - tc) 


(16) 


where 4>o is the pulse amplitude, Mqo is the free- 
st ream Mach number, which determines the width 
of the pulse, and tc is the nondimensional time at 
the center of the pulse. While a small pulse in 
roll angle is prescribed for the delta wing, the aero- 
dynamic transient is calculated. By using a transfer- 
function analysis, this aerodynamic transient is then 
used to obtain the rolling-moment coeflBcient in the 
frequency domain. In this case, a fast Fourier trans- 
form (FFT) of the rolling-moment coefficient is di- 
vided by an FFT of the pulsed rolling motion to ob- 
tain the value of Ci^. The transform assumes that 
the system is locally linear. Additional work not re- 
ported in this study has shown this assumption to be 
valid for the pulse amplitude of 1° used in the present 
study. 

Forced-Harmonic Analysis 

Because the pulse analysis is limited to small 
perturbations, the large-perturbation aerodynamic 
response characteristics of the rolling delta wing are 
determined using a forced- harmonic analysis. The 
forced-harmonic rolling motion can be expressed as 

(f) (t) — (j)o sin (ki) (17) 

where (po is now the roll amplitude, k is the reduced 
frequency of oscillation (based on half the wing root 
chord), and t is the nondimensional time. Since 
the linear techniques are no longer applicable, the 
concept of energy transferred to the system can 
be used in this analysis to determine the stability 
characteristics. A similar technique was applied 
experimentally in reference 8 for the analysis of wing- 
rock aerodynamics. During 1 cycle of harmonic 
motion, the total aerodynamic energy added to the 
system is 

AE= (fCi d(p (18) 
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where AE is a nondimensional energy and C/ is the 
rolling- moment coefficient. If AE > 0, then the aero- 
dynamic forces are. adding energy to the wing, which 
would have a destabilizing effect on the free-to-roll 
response. If AE < 0, then the aerodynamic forces 
are extracting energy from the wing, which would 
have a stabilizing effect on the free-to-roll response. 
Equation (18) indicates that for the rolling- moment 
coefficient versus roll-angle response, which traces a 
clockwise loop during 1 cycle of motion, the energy 
exchange is positive during the cycle. Similarly, for a 
counterclockwise loop, the energy exchange is nega- 
tive. If multiple loops are formed, then AE is a total 
of the individual loops. 

Free-to-Roll Analysis 


Roll Equation of Motion 

The equation of motion for a rolling delta wing 
can be expressed as 


— I ~ f^X^ ( 19 ) 

f 

where (p is the roll angle which is positive clockwise 
when viewed from aft, Ixx is the mass moment of 
inertia about the longitudinal axis, I is the aero- 
dynamic rolling moment, also positive clockwise, and 
is a structural damping term. (Dots over symbols 
indicate differentiation with respect to time.) To 
nondimensionalize equation (19), the angular rates 
are multiplied by the root chord of the delta wing c 
and divided by the free-stream speed of sound a^Q. 
The rolling-moment coefficient is defined as 


Cl 


I 

QooSc 


( 20 ) 


Time-Marching Solution 

The solution procedure for the time integration 
of equation (21) is based on a finite-difference repre- 
sentation of the time derivatives. The time deriva- 
tives are expressed in terms of second-order-accurate 
finite-difference approximations. After substituting 
these expressions into equation (21), the roll angle at 
time level n -h 1 can be expressed in terms of the roll 
angle at previous time levels as 


^ + (5 + 2C2Ai) r 

2C'2A? + 2 

(-4 - |C2AI) <^"-1 + ,^"-2 
^ |C'2A? + 2 


(23) 


The rolling moment at time level n'+l is esti- 

mated from a linear extrapolation of Ci at the pre- 
vious two time levels. This predicted value of Q 
is used to determine the roll angle at time 

level n + 1 . The flow field about the wing at this 
roll angle is then calculated, and the actual value of 
the rolling-moment coefficient is determined. The 
rolling-moment coefficient is then updated^ for use 
in the next time step. Because of the explicit 
time marching of the Euler code used in this study, 
the time steps required for numerical stability were 
small — approximately 6500 time steps per cycle of 
motion. Thus, it was not necessary to iterate between 
the roll angle calculation and the flow-field calcula- 
tion at each time step. For a free-to-roll calculation, 
steady-state initial conditions are specified for 
<tP, Cf\ a^nd Cf. At t = '0, an angular velocity 
perturbation is applied to the wing. 


where q^o is the free-stream dynamic pressure and 
S is the planform area. The nondimensional rolling 
equation of motion can then be written as 

/ = CiQ - (21) 


Active Roll Suppression 

Active roll suppression is achieved through the 
addition of an active rate-feedback control law to the 
time-marching solution procedure. A simple control 
law was chosen of the form 


where 


Cl 



C2 = 


^xC 

^oo^xx 


(22a) 

(22b) 


The structural damping term is added to simulate the 
damping that might be provided by a sting-balance 
bearing mount. This type of bearing mount was used 
in the low-speed wind-tunnel investigations of wing 
rock reported in references 8 to 11. 


6 = Ky<t>' (24) 

where Ky is the control gain and the values of <5 
are the left and right leading-edge flap deflection an- 
gles measured positive clockwise from the flap hinge 
lines. The control law is applied to the left and 
right flaps simultaneously, which results in an anti- 
symmetric configuration. The time-marching solu- 
tion procedure is the same as that described in the 
preceding section. However, after the roll angle at 
time level n-|-l is determined from equation (23), the 
flap deflection angle is determined from equation (24) 
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by using a second-order- accurate finite- difference ex- 
pression for the angular velocity (j)' . The deforming 
mesh algorithm is then applied, in addition to the 
rigid rotation, to move the mesh to its new posi- 
tion. As before, the flow field is calculated about 
the wing at its new position, and the rolling- moment 
coefficient is determined and updated for use in the 
next time step. The same initial conditions described 
in the preceding section are applied to begin the 
calculation. 

Results and Discussion 

Calculations were performed for a 75° swept delta 
wing at a free-stream Mach number of 1.2 and at 
a = 10°, 20°, and 30°. The wing has thickness and 
sharp leading edges as indicated in the partial view 
of the grid shown in figure 3. The thickness- to-span 
ratio at this cross section is 0.025, and the lower-edge 
bevel angle is 10°. The grid, which was generated 
by using :an advancing front method (ref. 45), has 
a total of 4226 nodes and 8299 elements. The grid 
was designed to be fine on the leeward side of the 
wing, where the dominant flow features are expected 
to occur. As discussed previously, the mesh is rotated 
as a rigid body for unsteady applications to conform 
to the instantaneous position of the main part of the 
wing. The mesh is deformed locally near the leading 
edges to conform to the instantaneous position of the 
flaps. As examples of mesh movement, partial views 
of the left leading-edge flap at a positive {6 = 10°) 
and a negative {6 = -10°) flap deflection angle are 
shown in figures 4(a) and 4(b), respectively, with 
the wing rotated through 10° of motion. The hinge 
point of the flap coincides with the inboard bevel 
edge on the lower surface; the flap length is therefore 
approximately 28 percent of the semispan. As shown 
in the figure, the mesh moves smoothly as the wing 
rolls and as the flaps are deflected. 

Steady and 'unsteady results, including the pulse, 
forced-harmonic, and free- to- roll calculations, are 
presented for a = 10°, 20°, and 30°. The rate- 
feedback control law is applied to the a = 30° case, 
since it is the only free-to-roll case to exhibit a wing- 
rock behavior. 

Steady-State Results 

Steady-state results were obtained to determine 
the basic character of the vortical flows and to pro- 
vide starting solutions for the unsteady cases. A com- 
parison of total pressure loss contours from these so- 
lutions (fig. 5) illustrates the effects of angle of attack. 
For the a = 10° case (fig. 5(a)), the contours indicate 
that the flow separates from each of the leading edges 
of the wing, which produces two small, widely spaced 


circular vortices. At a = 20° (fig. 5(b)), the contours 
indicate that the vortices are now larger than for the 
a = 10° case. The flow accelerating beneath the vor- 
tices at this increased angle of attack also produces 
two vertically oriented crossflow shock waves on the 
outboard portions of the wing. For the a — 30° case 
(fig, 5(c)), the contours indicate that the flow sep- 
arating from the leading edges produces two large, 
more closely spaced vortices. Also, as the flow ac- 
celerates beneath the vortices (fig. 5(c)), vertically 
oriented crossflow shock waves are formed on the out- 
board portions of the wing. Weaker shock waves are 
formed on the top of each vortex. These vertically 
oriented shocks are located above the core of the vor- 
tex. A weaker horizontal shock wave is also present 
between the vortices. 

Pulse Transfer-Function Results 

The pulse transfer-function analysis was per- 
formed to determine the small amplitude stability 
and response characteristics of the wing. A sta- 
bility analysis was derived by first recalling the 
nondimensional rolling equation of motion given by 

cP" = CiCi - C2<j>' (21) 

For simplicity, the structural damping term C 2 is set 
equal to zero, which results in 

f = CiCi (25) 

Assuming that the rolling-moment coefficient can 
be written as the product of the rolling- moment- 
coefficient transfer function and the roll angle 0, 
then 

Cl = Ci^<p (26) 

Then, for simple harmonic motion, 

Ci = [Re{Ci^) + ilm{Ci^)](t> (27) 

where Re(C/^) and Im(C/^) represent the real and 

imaginary parts of the first harmonic component of 
Ci^. In this case, the real part of the rolling-moment- 
coefficient transfer function represents an aero- 
dynamic stiffness, and the imaginary part represents 
an aerodynamic damping. Therefore, for this sim- 
ple one-degree-of- freedom case, the sign of Im(Q^) 
determines the stability in roll of the wing for small 
perturbations. In other words, a negative value of 
Im(C'/^) indicates a positive aerodynamic damping, 
which would cause a free-to-roll wing to be stable; 
a positive Im(Q^) indicates a negative aerodynamic 
damping, which would cause a free-to-roll wing to be 
unstable. 
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The pulse transfer-function analysis is used to de- 
termine the force-coefficient transfer function, and 
therefore the stability of the wing, for a wide range of 
reduced frequency k. The nondimensional time step 
used in these analyses was At = 0.004, which results 
in a frequency resolution of A A: — 0.1. A comparison 
of results from the pulse transfer-function analysis 
(fig. 6) indicates the effects of angle of attack. A max- 
imum pulse amplitude of 1° was used at each angle 
of attack. However, at o: = 20°, an additional pulse 
amplitude of —1° was also considered for reasons ex- 
plained below. At a = 10° (fig. 6(a)), Im(Q^) is 
negative for all values of reduced frequency, which is 
indicative of stability in roll for small perturbations. 
At a = 20° (fig. 6(b)), Im(C'/^) is also negative over 
all values of reduced frequency for both pulse angles. 
However, the absolute value of Im(C';^) for each value 
of k is smaller for a = 20° results than for a = 10° 
results; thus, the aerodynamic damping should be 
lower at a = 20°. In contrast, at o; — 30° (fig. 6(c)), 
the imaginary part is positive for k < 0.5, which is 
indicative of instability in roll. Also, the Re(Q^) 
is negative for a = 30° in this range of reduced fre- 
quency, which corresponds to a positive aerodynamic 
stiffness. The roll response will therefore oscillate 
(with increasing amplitude) rather than give rise to 
a static instability known as wing drop. 

The accuracy of the pulse analysis is verified by 
harmonic analyses performed at five values of re- 
duced frequency: k = 0.0, 0.25, 0.50, 0.75, and 1.0. 
In these analyses, the wing was oscillated harmoni- 
cally in roll with an amplitude of 1° for three cycles 
of motion; the rolling-moment coefficient was deter- 
mined from the last cycle. The nondimensional time 
step used in all cases was 0.00262, which corresponds 
to 2600 steps per cycle for the k = 0.25 cases and 
650 steps per cycle for the A: = 1.0 cases. The re- 
sults of the harmonic analyses are compared with the 
rolling-moment-coefficient transfer functions shown 
in figure 6. The agreement between pulse and har- 
monic analyses is very good for a — 10° and 30°. For 
a = 20°, the harmonic results lie between the pulse 
results for 4>o = 1° and (j)o = — . Because of the ini- 

tial symmetry of the steady-state configuration, the 
pulse results should be independent of the sign of the 
pulse amplitude. However, flow- field asymmetries as- 
sociated with the asymmetric mesh cause some bias 
in the pulse-analysis results at this angle of attack 
that are not present in the harmonic analyses, be- 
cause the harmonic analyses involve oscillations be- 
tween (j) = ±1°. This symmetric motion effectively 
compensates for the flow-field asymmetries. 


Forced-Harmonic Results 

Because the pulse transfer-function analysis is 
limited to small perturbations, the large-perturbation 
aerodynamic characteristics of the delta wing were 
investigated with forced-harmonic oscillations. Mo- 
tions at a reduced frequency of A; = 0.25 were chosen 
for this analysis. This value lies at the midpoint of 
the range of reduced frequency that is identified by 
the pulse analysis as being an unstable condition for 
the free-to-roll wing at o = 30°. Three amplitudes of 
motion — = 5°, 15°, and 35° — were considered at 
a = 10°, 20°, and 30°. The nondimensional time step 
used for all cases was 0.00262, which corresponds to 
4000 steps per cycle of harmonic motion. A compari- 
son of rolling-moment coefficient versus roll angle for 
each of these cases is shown in figure 7 to illustrate 
the eflPects of roll amplitude and angle of attack. For 
a = 10° (fig. 7(a)), the results indicate a counter- 
clockwise loop for each roll amplitude, which would 
produce a convergent (stable) response if the wing 
were free to roll. This prediction of a stable response 
at the smallest roll amplitude is consistent with the 
pulse transfer-function results in figure 6(a). Also, 
as the roll amplitude is increased from 5° to 15°, the 
maximum rolling-moment coefficient increases lin- 
early. (Note the change in scaling of the vertical 
and horizontal axes.) However, as the roll ampli- 
tude is further increased to 35°, some nonlinear aero- 
dynamic characteristics are exhibited in the “pinch- 
ing” of the loop at the extreme roll angles; however, 
the free-to-roll response is still predicted to be stable. 
Similarly, for a = 20° (fig. 7(b)), the results show a 
counterclockwise loop (A£ < 0) for each roll ampli- 
tude, which also would produce a stable response if 
the wing were free to roll. This prediction of a stable 
response at the smallest roll amplitude, which would 
be more lightly damped than for a = 10°, is also con- 
sistent with the pulse transfer-function results of fig- 
ure 6(b). For a = 30° (fig. 7(c)), the results indicate 
clockwise loops for the 5° and 15° roll amplitudes, 
which would produce a divergent (unstable) response 
if the wing were free to roll. This prediction of an 
unstable, free-to-roll response at the smaller roll am- 
plitudes is consistent with the pulse transfer-function 
results of figure 6(c). For ~ 35°, counterclockwise 
loops have formed at the extreme roll angles which, 
consequently, would have a stabilizing effect on the 
free-to-roll response. The formation of these stabiliz- 
ing loops was not, of course, predicted by the pulse 
analysis. In contrast to the a = 10° and 20° cases, 
the nonlinear aerodynamic effects at the larger roll 
amplitudes for a = 30° result in a change in the 
stability characteristics of the wing. 
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The total aerodynamic energy exchange during 
a cycle of motion, as described in equation (18), 
was calculated for each roll amplitude. These non- 
dimensional energy values were then normalized by 
a factor of 0 oA^ 3 O, where (j>o is the roll amplitude 
during the cycle (in degrees) and is the energy 

exchange value for a = 30° and (f)o = 1^^- The 
normalized energy exchange during the cycle can 
thus be expressed as 



Plots of versus (f)o at each angle of attack are 
shown in figure 8. Additional results for the roll am- 
plitudes of 1°, 30°, and 45° are included in figure 8. 
If the aerodynamic effects were linear ov^ the range 
of roll amplitudes, then the curve of AE ve^us (po 
for a = 30° should be a horizontal line at A£^ = 1. 
Similarly, the curves for a = 10° and 20° should be 
horizontal lines at their respective normalized energy 
exchange values for (po — 1° if the system was lin- 
ear. For a = 10° and 20°, the near- horizontal lines 
shown in figure 8 indicate that the aerodynamic ef- 
fects are nearly linear with respect to roll amplitude. 
For a = 30° , nonlinear effects are indicated in figure 8 
at the larger roll amplitudes, which would result in a 
change in the stability characteristics of the free-to- 
roll wing. A neutr^ally stable condition (AE = 0) or 
limit cycle oscillation is predicted at approximately 
(p^ = 36° for this reduced frequency. 

The unsteady vortex dynamics during a harmonic 
cycle can be illustrated by the changes in the cross- 
flow contours during the cycle. Crossflow total- 
pressure-loss contours from the a — 30° forced- 
harmonic cases at (po = 5° and 35° are shown in 
figure 9. The instantaneous crossflow contours are 
shown at four points in time that correspond to the 
0° (1), 90° (2), 180° (3), and 270° (4) cycle posi- 
tions. For (po = 5° (fig. 9(a)), the contours show only 
a slight variation in the vortex strength and location 
during the cycle. However, for (po = 35° (fig. 9(b)), 
the contours indicate that the vortex strength and 
location change significantly during the cycle. Fig- 
ure 9(b) shows that as the left leading edge moves 
through zero roll angle (position 1) and continues 
to the maximum position (position 2), the left vor- 
tex weakens and lifts off the wing, while the right 
vortex strengthens and moves inboard. Similarly, as 
the right leading edge moves through zero roll angle 
(position 3) to the maximum position (position 4), 
the right vortex weakens and lifts off the wing, while 
the left vortex strengthens, reattaches, and moves 
inboard. This vortex lift-off and reattachment is be- 
lieved to be the source of the nonlinear variation 


of the rolling-moment coefficient shown in figures 7 
and 8. 

Free-to-Roll Results 

The free-to-roll results were obtained for the flow 
conditions and structural and inertial parameter val- 
ues listed in table 1. The structural and inertial 


Table 1. Summary of Structural Parameter Values and 
Flow Conditions for Free-to-Roll Calculation 


c, m . . . . 
Ixx, kg-m^ . 
//x, kg-m^/sec 
Poo, kg/m^ . 
aoo, m/sec . 


, . . . 0.282 

0.1776 X 10’^ 

, 0 

. . . . 0.526 

312 


properties used in these calculations are loosely based 
on the characteristics of the models used in the ex- 
perimental study of wing rock in reference 9. The 
initial nondimensional angular velocity imposed on 
the wing was 0.003, and the nondimensional time 
step was 0.004. The resulting roll- angle response 
for the a = 10° case (fig. 10(a)) indicates that, af- 
ter the initial perturbation, the oscillatory response 
converges to its initial steady- state value. This sta- 
ble free-to-roll response is consistent with the pulse 
and forced-harmonic results presented in figures 6(a) 
and 7(a). Similarly, the resulting roll- angle response 
for a = 20° (fig. 10(b)) also shows a stable, con- 
verging response that is consistent with the pulse 
and forced-harmonic results of figures 6(b) and 7(b). 
Also, as predicted by the pulse and harmonic analy- 
ses, the response at a = 20° is more lightly damped 
than the response at a = 10°. The roll-angle re- 
sponse for a = 30° (fig. 10(c)) indicates that the 
oscillatory response initially diverges for small val- 
ues of roll angle, which is consistent with the small 
amplitude pulse and harmonic results of figures 6(c) 
and 7(c). As the roll angle increases to around 35°, 
the rate of divergence decreases because of the sta- 
bilizing aerodynamics (counterclockwise loops in the 
rolling-moment coefficient at the extreme roll angles) 
shown in figure 7(c). Finally, the response reaches 
a maximum amplitude of motion at approximately 
(p — 38°; this response corresponds to a limit cycle 
oscillation. The reduced frequency of the limit cycle 
is 0.103. 

In the low-speed experimental investigation of 
wing rock conducted by Arena and Nelson (ref. 11), 
the wing-rock time history of an 80° swept delta 
wing at 30° angle of attack exhibited a symmet- 
ric growth of the maximum roll angle which, at the 
larger roll amplitudes, transitioned to a limit cycle 
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oscillation at fc = 0.125. During the wing-rock mo- 
tion, no vortex burst was observed above the config- 
uration. For the high-speed calculations performed 
in this study, the wing-rock time history exhibited a 
similar symmetric growth and transition. However, 
the structured-grid, conical Euler calculation from 
reference 33 for an 80° swept delta wing at a free- 
stream Mach number of 1.2 and 35° angle of attack 
indicated an antisymmetric growth of the maximum 
roll amplitude that transitioned to a limit cycle os- 
cillation about a mean roll angle of —5°. 

Similar to the forced- harmonic results shown in 
figure 9, the unsteady vortex dynamics during the 
wing-rock cycle are illustrated in figure 11 by the 
changes in the crossflow total-pressure-loss contours 
during the cycle. The instantaneous crossflow con- 
tours are shown at four points in time during the 
wing-rock cycle, these points correspond to the 
0° (1), 90° (2), 180° (3), and 270° (4) cycle posi- 
tions. Figure 11 shows that, similar to the results 
shown in figure 9(b), as the left leading edge moves 
through zer6 roll angle (position 1) and continues to 
the maximum position (position 2), the left vortex 
weakens and lifts off the wing and the right vortex 
strengthens and moves inboard. As the right leading 
edge moves up to its maximum position (position 4), 
the right vortex weakens and lifts off the wing and 
the left vortex strengthens, reattaches, and moves 
inboard. The similarities between the results of fig- 
ure 9(b) and figure 11 lead to the conjecture that 
the vortex lift-off and reattachment are the sources 
of the change in aerodynamic damping that stabi- 
lizes the wing response; thus, a limit cycle oscillation 
is produced. However, the details of the fluid mech- 
anisms that produce the wing rock are still under 
investigation. 

It is important to address two of the major lim- 
itations of the conical Euler equations. The first of 
these is that the conical Euler equations cannot pre- 
dict the formation of secondary vortices. However, 
it is shown in reference 24 that for steady flow, the 
Euler equations accurately model the primary vortex 
for a sharp-edge delta wing. The second limitation 
of the conical Euler equations is that these equations 
cannot predict the time lag in the radial directions. It 
was previously noted that the unsteady conical Euler 
equations assume instantaneous propagation of dis- 
turbances in the radial directions. However, as the 
free-stream Mach number increases and the reduced 
frequency of oscillation decreases, this approximation 
improves. Although the cases presented in this work 
are not for extremely high Mach numbers, the re- 
duced frequency of oscillation is low. For example, at 
a reduced frequency of 0.103, the phase shift from the 


apex to the trailing edge is approximately 7°. The 
effects of this phase lag cannot at this time be accu- 
rately quantified. However, if an analogy is made be- 
tween the conical calculations and the computational 
work done on two-dimensional airfoils, then the con- 
ical results can be thought of in terms of crossflow 
sectional properties of the delta wing. These coni- 
cal studies can then predict qualitative information 
about the sectional properties of the delta wing at dif- 
ferent flow conditions, which can give insight and di- 
rection to subsequent three-dimensional calculations. 
Also, as in this study, the conical equations can act as 
an efficient test-bed for developing analysis methods 
that can be directly extended to three-dimensional 
calculations. 

Active Roll-Suppression Results 

An active rate-feedback control law was imple- 
mented in an attempt to suppress the wing-rock mo- 
tion. To determine an appropriate value for the gain, 
a stability analysis was derived by again using the 
nondimensional rolling equation of motion given by 


= CiQ (25) 

Assuming that the rolling-moment coefficient can be 
written as the superposition of the rolling-moment- 
coefficient transfer functions for (j) and 5, then 


Ci = Ci^<f> + Ci^6 (29) 

Substituting the control law from equation (24) into 
equation (29) gives 

Q = + (30) 

Then, for simple harmonic motion. 

Cl = {[Re(Q^) - fcMooR. Im(Q,)] ] 

> (31) 

+ i[lm{Ci^) + fcMooR. Re(Q,)]} J 

As before, the first term on the right-hand side 
of equation (31) represents an aerodynamic stiff- 
ness, and the second term represents an aero- 
dynamic damping. Therefore, stabilizing the motion 

of the wing requires that the aerodynamic damping 
be positive, or that 


Im(Q^) + kMopKy Re(Q^) < 0 
Solving for the gain yields 


Ky > 


-1 

fcMoo Re(Q^) 


(32) 


(33) 
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A pulse analysis was performed to determine the 
rolling-moment-coefficient transfer function for 6 at 
a = 30° (Pulse amplitude = 1°). The transfer func- 
tion from this analysis (fig. 12) indicates that for 
values of reduced frequency less than 0.5, the real 
part of is negative. Considering equation (33) 
and recalling from figure 6(c) that the values of 
Im(Q^) for k < 0.5 are positive, the value of Ky 
must be positive to suppress wing rock. The value 
of Ky actually needs to be greater than that de- 
termined by evaluating the right-hand side of equa- 
tion (33) to stabilize the wing, since the preceding 
analysis assumes simple harmonic motion. For the 
flow conditions considered herein, the value for the 
gain that produces a neutrally stable (or simple har- 
monic) response is Ky =0.35. The free- to- roll anal- 
ysis with active rate-feedback control was performed 
for Ky = 0.25, 0,40, and 0.50. The time histories 
of the wing motion are shown in figure 13. As ex- 
pected, the time history for Ky = 0.50 (fig. 13(a)) 
indicates a damped response. Similarly, the response 
for Ky = 0.40 (fig. 13(b)) is also damped, although 
at a smaller rate than for Ky = 0.50. The response of 
the wing for Ky = 0.25 (fig. 13(c)) indicates that the 
response is no longer damped. However, a compar- 
ison with the results of figure 10(c) {Ky = 0) shows 
that the active rate-feedback control reduces the 
growth rate of the response. The maximum flap de- 
flection commanded by the control law during these 
free-response calculations was in the range of 1° to 2°. 

Concluding Remarks 

A conical Euler analysis method was developed 
to study unsteady, vortex-dominated flows about 
rolling, highly swept delta wings undergoing either 
forced motions or free-to-roll motions that include 
active roll suppression. The flow solver of the code 
involves a multistage, Runge-Kutta time-stepping 
scheme that uses a cell- centered, finite-volume, spa- 
tial discretization of the Euler equations on an un- 
structured grid of triangles. The code allows for the 
additional analysis of the free-to-roll case by simulta- 
neously integrating in time the rigid-body equation 
of motion with the governing flow equations. Results 
are presented for a delta wing with a 75° swept, sharp 
leading edge at a free-stream Mach number of 1.2 
and at 10°, 20°, and 30° angle of attack a. At the 
lower angles of attack (10° and 20°), forced- harmonic 
analyses indicate that the rolling-moment coefficients 
provide a positive damping, which is verified by free- 
to-roll calculations. In contrast, at the higher angle 
of attack (30°), a forced-harmonic analysis indicates 
that the rolling-moment coefficient provides negative 
damping at the small roll amplitudes. A free-to-roll 


calculation for this case produces an initially diver- 
gent response, but as the amplitude of motion grows 
with time, the response transitions to a wing-rock 
type of limit cycle oscillation, which is characteris- 
tic of highly swept delta wings. The wing-rock mo- 
tion was subsequently suppressed, however, by the 
use of an active rate-feedback control law and anti- 
symmetrically deflected leading-edge flaps. The 
methodology developed is directly extendable to 
three-dimensional calculations. 
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Figure 1. Example of a central-difference type of differencing approach. 
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Figure 2. Difference stencil for harmonic operator. 



Figure 3. Partial view of unstructured grid about 75° swept delta wing. 
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(a) 8 = 10°. 



(b) <5 = -10°. 

Figure 4. Partial view of deforming mesh about deflected leading-edge flap. 
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0.04 5 P{ S 0.56 
Apj = 0.04 



(b) a = 20°. 



0.04 < P( < 0.76 
Apj = 0.04 


(c) a = 30°. 

Figure 5. Effects of angle of attack on steady-state total-pressure-loss contours for a 75° swept delta wing at 
Moo = 1-2. 
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(c) a = 30^. 

Figure 6. Effects of angle of attack on rolling-moment-coefficient transfer function versus reduced frequenc} 
for a 75° swept delta wing at Moo = 1-2. 








(a) 4>o = b°. 

Figure 9. Total-pressure-loss contours during a cycle of harmonic rolling motion for a 75° swept delta wing at 
Moo = 1.2, a = 30°, and k = 0.25. 
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(b) 4>o = 35°. 
Figure 9. Concluded. 
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